/*
 * ReactionSystem.cpp
 *
 *  Created on: 23 Jul 2011
 *      Author: Allan
 */

#include "ReactionSystem.h"

// C++ includes
#include <iomanip>
#include <set>

// Boost includes
#include <boost/foreach.hpp>

// Eigen includes
#include <Eigen/Core>
using namespace Eigen;

ReactionSystem::ReactionSystem()
{}

ReactionSystem::ReactionSystem(const vector<Reaction>& reactions) :
reactions(reactions)
{}

void ReactionSystem::ArrangeReactions(const vector<string>& reactants)
{
	auto index = [&](const string& s) { return std::find(reactants.begin(), reactants.end(), s) - reactants.begin(); };
	
	auto compare = [&](const Reaction& r1, const Reaction& r2) { return index(r1.GetReactant().name) < index(r2.GetReactant().name); };
	
	std::sort(reactions.begin(), reactions.end(), compare);
}

///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
/// Basis Swap Implementation todo finish this implementation
///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//void ReactionSystem::RestructureReactions(const vector<string>& newPrimarySpecies, const vector<string>& newSecondarySpecies)
//{
//	uint Nj = products.size();
//	uint Ni = reactants.size();
//	
//	MatrixXd M = MatrixXd::Zero(Nj+Ni, Nj);
//	
//	M.block( 0, 0, Nj, Nj) = MatrixXd::Identity(Nj, Nj);
//	M.block(Nj, 0, Ni, Nj) = vij;
//	
//	// The primary-swap matrix (theta)
//	MatrixXd theta = MatrixXd::Zero(Nj, Nj);
//	
//	// Loop over all the species in (newPrimarySpecies) and set the primary-swap matrix (theta) 
//	for(uint row = 0; row < Nj; ++row)
//	{
//		// Determine the local indexes of the current primary species in the set of the actual primary and secondary species, (primarySpecies) and (secondarySpecies) respectively (note that only one of the indexes below will be within range) 
//		const unsigned j = std::find(products.begin(), products.end(), newPrimarySpecies[row]) - products.begin();
//		const unsigned i = std::find(reactants.begin(), reactants.end(), newPrimarySpecies[row]) - reactants.begin();
//		
//		// Check if the current new primary species (primary) is in the set of primary species (primarySpecies) and set the row of (theta) accordingly 
//		if(j < Nj) theta(row, j) = 1.0;
//		
//		// Check if the current new primary species (primary) is in the set of secondary species (secondarySpecies) and set the row of (theta) accordingly 
//		if(i < Ni) theta.row(row) = vij.row(i);
//	}
//	
//	// Transform the matrix M to reflect the changes in the set of primary species (M' = M * inv(theta))
//	MatrixXd thetaInv = FullPivLU<MatrixXd>(theta).inverse();
//	
//	M *= thetaInv;
//	
//	// Update the stoichiometric matrix vij
//	for(uint row = 0; row < Ni; ++row)
//	{
//		// Determine the local indexes of the current primary species in the set of the actual primary and secondary species, (primarySpecies) and (secondarySpecies) respectively (note that only one of the indexes below will be within range) 
//		const unsigned j = std::find(products.begin(), products.end(), newSecondarySpecies[row]) - products.begin();
//		const unsigned i = std::find(reactants.begin(), reactants.end(), newSecondarySpecies[row]) - reactants.begin();
//		
//		vij.row(row) = (j < Nj) ? M.row(j) : M.row(i + Nj);
//	}
//	
//	// Create the new reactions
//	vector<Reaction> newReactions(Ni);
//	
//	for(uint i = 0; i < Ni; ++i)
//	{
//		vector<string> products;
//		vector<double> stoichiometries;
//		
//		const double tol = 1.0E-6;
//		
//		for(uint j = 0; j < Nj; ++j)
//		{
//			if(std::abs(vij(i, j)) < tol)
//			{
//				products.push_back(newPrimarySpecies[j]);
//				stoichiometries.push_back(vij(i, j));
//			}
//		}
//		
//		newReactions[i].SetEquation(newSecondarySpecies[i], products, stoichiometries);
//	}
//	
//	
//}

void ReactionSystem::SetIndexes(const SpeciesMap& speciesMap)
{
	// Loop over all the reactions and set their index related data
	BOOST_FOREACH(Reaction& reaction, reactions)
		reaction.SetIndexes(speciesMap);
}

const vector<Reaction>& ReactionSystem::GetReactions() const
{
	return reactions;
}

const vector<string> ReactionSystem::Reactants() const
{
	// The reactants of the reactions
	vector<string> reactants;
	
	// Loop over all the reactions
	BOOST_FOREACH(const Reaction& reaction, reactions)
		reactants.push_back(reaction.GetReactant().name);
	
	return reactants;
}

const vector<string> ReactionSystem::Products() const
{
	// The unique products of the reactions
	set<string> products;
	
	// Loop over all the reactions and insert their products into the set of products
	BOOST_FOREACH(const Reaction& reaction, reactions)
		BOOST_FOREACH(const Reaction::Species& product, reaction.GetProducts())
			products.insert(product.name);

	return vector<string>(products.begin(), products.end());
}

const MatrixXd ReactionSystem::StoichiometricMatrix() const
{
	return StoichiometricMatrix(Products());
}

const MatrixXd ReactionSystem::StoichiometricMatrix(const vector<string>& products) const
{
	// The number of reactions (Nr) and products (Nj)
	const uint Nr = reactions.size();
	const uint Nj = products.size();

	// The stoichiometric matrix for the set of reactions
	MatrixXd vrj = MatrixXd::Zero(Nr, Nj);

	// Loop over all the reactions
	for(uint r = 0; r < Nr; ++r)
	{
		// Loop over all the products in the current reaction
		BOOST_FOREACH(const Reaction::Species& product, reactions[r].GetProducts())
		{
			// Determine the index of the the current (product) in the vector of (products)
			const uint j = std::find(products.begin(), products.end(), product.name) - products.begin();
			
			// Set the entry (r,j) of the stoichiometric matrix
			vrj(r, j) = product.stoichiometry; 
		} 
	}
	
	return vrj;
}

const VectorXd ReactionSystem::EquilibriumConstants(double T, double P) const
{
	VectorXd Ki(reactions.size());
	
	for(uint i = 0; i < (uint)Ki.rows(); ++i) Ki[i] = reactions[i].EquilibriumConstant(T, P);
	
	return Ki;
}

const VectorXd ReactionSystem::ReactionQuotients(const VectorXd& a) const
{
	VectorXd Qi(reactions.size());
	
	for(uint i = 0; i < (uint)Qi.rows(); ++i) Qi[i] = reactions[i].ReactionQuotient(a);
	
	return Qi;
}

const VectorXd ReactionSystem::MassRates(double T, double P, const VectorXd& n, const VectorXd& a) const
{
	VectorXd Ri(reactions.size());
	
	for(uint i = 0; i < (uint)Ri.rows(); ++i) Ri[i] = reactions[i].MassRate(T, P, n, a);
	
	return Ri;
}

void ReactionSystem::PrintState(double T, double P, const VectorXd& n, const VectorXd& a) const
{
	cout << setw(25) << left << "Reactants";
	cout << setw(25) << left << "log(Qr)";
	cout << setw(25) << left << "log(Kr)";
	cout << setw(25) << left << "log(Kr/Qr)" << endl << endl;
	
	BOOST_FOREACH(const Reaction& r, reactions)
	{
		const double Qr = r.ReactionQuotient(a);
		const double Kr = r.EquilibriumConstant(T, P);
		
		cout << setw(25) << left << r.GetReactant().name;
		cout << setw(25) << left << std::log10(Qr);
		cout << setw(25) << left << std::log10(Kr);
		cout << setw(25) << left << std::log10(Kr/Qr) << endl;
	}
}

ostream& operator<<(ostream& out, const ReactionSystem& reactionSystem)
{
	out << "Reactions: " << endl;
	
	BOOST_FOREACH(const Reaction& reaction, reactionSystem.GetReactions())
		out << '\t' << reaction << endl;
	
	out << endl;
	
	out << "Reactants: " << endl << '\t';
		out << reactionSystem.Reactants() << endl << endl;
		
	out << "Products: " << endl << '\t';
		out << reactionSystem.Products() << endl << endl;
		
	out << "Stoichiometric Matrix: " << endl;
		out << reactionSystem.StoichiometricMatrix(reactionSystem.Products()) << endl;
		
	return out;
}

